#include "config.h"

#include "lstatfit.h"

#define GAP	0.00001

ecm_double
weibull_goodness(msl_fit_t *fit)
{
	ecm_double	*xvalues = fit->xvalues;
	ecm_double	*yvalues = fit->yvalues;
	ecm_uint32	i;
	double	k, lambda;
	ecm_double	sq_diffs = 0;

	k = gsl_vector_get(fit->s->x, 0);
	lambda = gsl_vector_get(fit->s->x, 1);

	for (i = 0; i < fit->n; i++) {
		double	t = (xvalues[i] + GAP) / (fit->xrange + 2 * GAP) * fit->xrange;
		ecm_double	prob;

		prob = 1 - exp(-pow(t / lambda, k));
		sq_diffs += pow(yvalues[i] - prob, 2);
	}

	return sq_diffs;
}

void
weibull_dist_report(FILE *fp, msl_fit_t *fit, ecm_uint32 xdiff)
{
	ecm_double	*xvalues = fit->xvalues;
	ecm_uint32	i;
	double	k, lambda;

	k = gsl_vector_get(fit->s->x, 0);
	lambda = gsl_vector_get(fit->s->x, 1);

	ecm_fprintf(fp, "# k:%lf lambda:%lf\n", k, lambda);

	for (i = 0; i < fit->n; i++) {
		double	t = (xvalues[i] + GAP) / (fit->xrange + 2 * GAP) * fit->xrange;
		ecm_double	prob;

		prob = 1 - exp(-pow(t / lambda, k));
		ecm_fprintf(fp, "%lf %lf\n", xvalues[i] + xdiff, prob);
	}
}

int
weibull_dist_f(const gsl_vector *x, void *data, gsl_vector *f)
{
	msl_fit_t	*fit = (msl_fit_t *)data;
	ecm_double	*xvalues = fit->xvalues;
	ecm_double	*yvalues = fit->yvalues;
     
	double	k = gsl_vector_get(x, 0);
	double	lambda = gsl_vector_get(x, 1);

	size_t	i;
     
	for (i = 0; i < fit->n; i++)	{
		/* Model Yi = 1 - exp(-(t/lambda)^k) */
		double	t = (xvalues[i] + GAP) / (fit->xrange + 2 * GAP) * fit->xrange;
		double	Yi = 1 - exp(-pow(t / lambda, k));

		gsl_vector_set(f, i, Yi - yvalues[i]);
	}
     
	return GSL_SUCCESS;
}

int
weibull_dist_df(const gsl_vector *x, void *data, gsl_matrix *J)
{
	msl_fit_t	*fit = (msl_fit_t *)data;
	ecm_double	*xvalues = fit->xvalues;
	ecm_uint32	i;
	double	k = gsl_vector_get(x, 0);
	double	lambda = gsl_vector_get(x, 1);

	for (i = 0; i < fit->n; i++) {
		double	t = (xvalues[i] + GAP) / (fit->xrange + 2 * GAP) * fit->xrange;
		double	dk;
		double	dlambda;

		dk = exp(-pow(t / lambda, k)) * pow(t / lambda, k) * log(t / lambda);
		dlambda = -exp(-pow(t / lambda, k)) * k * pow(t, k) / pow(lambda, k + 1);
		
		gsl_matrix_set(J, i, 0, dk);
		gsl_matrix_set(J, i, 1, dlambda);
	}
	return GSL_SUCCESS;
}

